Jedną z możliwych metod stabilizacji równań Stokesa jest metoda DG (Discontinuous Galerkin). W metodzie tej funkcje bazowe używane do aproksymacji i testowania rozwiązania "rozcinane" są nożyczkami na granicy poszczególnych elementów. Tak więc jeśli chcemy używać wielomianów B-spline, możemy je dalej używać, ale nad każdym elementem zostaną one rozcięte na lokalne kawałki.
W rozdziale tym opieramy się na sformułowaniu DG równań Stokesa na podstawie książki [1]
Używamy sformułowania wyprowadzonego w poprzednim rozdziale
\( a({\bf u},{\bf v}) + b(p,{\bf v }) = f({\bf v}) \quad \forall {\bf v} \)
\( b(q,{\bf u } ) = 0 \quad \forall q \)
\( a({\bf u}, {\bf v}) =\int_{\Omega} \nabla {\bf u} : \nabla {\bf v }dxdy = \int_{\Omega} \frac{\partial u_1}{\partial x}\frac{\partial v_1}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_1}{\partial y}\frac{\partial v_1}{\partial y}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial x}\frac{\partial v_2}{\partial x}dxdy + \int_{\Omega} \frac{\partial u_2}{\partial y }\frac{\partial v_2}{\partial y }dxdy \)
\( b(p, {\bf v}) = \int_{\Omega} p div {\bf v} dxdy= \int_{\Omega}p \frac{\partial v_1}{\partial x} dxdy+ \int_{\Omega} p \frac{\partial v_2}{\partial y } dxdy \)
\( f({\bf v} )= \int_{\Omega} \frac{\partial h(x,y)}{\partial x}\frac{\partial v_1(x,y)}{\partial x } dxdy +\int_{\Omega} \frac{\partial h(x,y)}{\partial y}\frac{\partial v_2(x,y)}{\partial y} dxdy +\int_{\partial \Omega} \frac{\partial h(x,y)}{\partial n } v_1(x,y) dS \)
W przypadku sformułowania dla metody DG konieczne jest wprowadzenie dodatkowych członów opisujących skoki funkcji bazowych
\( [\psi(x)]=\psi|_{T_1}(x)-\psi|_{T_2}(x) \) gdzie \( T_1 \) i \( T_2 \) oznacza dwa sąsiednie elementy, i obliczamy tutaj różnce z wartości "pociętej" funkcji B-spline oznaczonej \( \psi \) na elemencie z jednej strony i z drugiej strony.
Dodatkowo wprowadza się średnią z wartości rozciętych funkcji bazowych na krawędziach dwóch sąsiednich elementów \( \{\psi\} = \frac{\psi|_{T_1}(x)+\psi|_{T_2}(x)}{2} \).
Nasz funkcjonał \( a({\bf u },{\bf v }) \) należy zmodyfikować dla metody DG w następujący sposób
\( a_{DG}({\bf u }, {\bf v }) = \sum_{K \in \Omega_h } \int_K \left( \frac{\partial u_1 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_1 }{\partial y } \frac{\partial v_1 }{\partial y } + \frac{\partial u_2 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_2 }{\partial y } \frac{\partial v_1 }{\partial y } \right) \)
\( + \sum_{F\in F_h} \int_F \frac{\eta}{h_F} \left( [\frac{\partial u_1}{\partial x} ] [\frac{\partial v_1}{\partial x}] + [\frac{\partial u_1}{\partial y}] [\frac{\partial v_1}{\partial y }] + [\frac{\partial u_2}{\partial x} ] [\frac{\partial v_2}{\partial x}] + [\frac{\partial u_2}{\partial y}] [\frac{\partial v_2}{\partial y }]\right) \)
\( -\sum_{F\in F_h} \int_F \left( \{\frac{\partial u_1}{\partial n_1} \} [\frac{\partial v_1}{\partial x}] + \{ \frac{\partial u_1}{\partial n_2} \} [\frac{\partial v_1}{\partial y } ] + \{\frac{\partial u_2}{\partial n_1} \} [\frac{\partial v_2}{\partial x}] + \{ \frac{\partial u_2}{\partial n_2} \} [\frac{\partial v_2}{\partial y } ]\right) \)
\( - \sum_{F\in F_h} \int_F \left( [\frac{\partial u_1}{\partial x} ] \{ \frac{\partial v_1}{\partial n_1} \} + [\frac{\partial u_1}{\partial y} ] \{\frac{\partial v_1}{\partial n_2 } \} + [\frac{\partial u_2}{\partial x} ] \{ \frac{\partial v_2}{\partial n_1} \} + [\frac{\partial u_2}{\partial y} ] \{\frac{\partial v_2}{\partial n_2 } \} \right) \)
W powyższych wzorach, całki liczone są element po elemencie. Tak więc \( \Omega_h \) oznacza zbiór wszystkich elementów z patcha na którym rozpięte są nasze funkcje B-spline, \( K \in \Omega_h \) oznacza pojedynczy element wybrany z rodziny elementów (tak więc sumujemy nasze całki elementowe po wszystkich elementach z patcha), \( F_h \) oznacza zbiór wszystkich krawędzi pomiędzy elementami patcha, \( F\in F_h \) oznacza pojedynczą krawędź wybraną z tej rodziny (tak więc sumujemy nasze całki krawędziowe po wszystkich krawędziach z patcha), \( h_F \) oznacza dugość krawędzi \( F \), natomiast \( \eta \) oznacza parametr którego wartość należy dobrać eksperymentalnie.
Niestety wadą tej metody stabilizacji jest fakt, iż musimy uruchomić kilka obliczeń i znaleźć taką wartość parametru \( \eta \) dla której symulacja da poprawne wyniki. Zazwyczaj eksperymentuje się z wielkościami różniącymi się o rząd wielkości np. 0.1, 1, 10, 100, 1000, 10000 etc.
Symbol \( {\bf n }=(n_1,n_2) \) oznacza wersor prostopadły do krawędzi \( F \). W naszym przypadku kiedy rozważamy funkcje bazowe B-spline rozpięte na regularnym patchu elementów w 2D, przyjmujemy \( {\bf n }=(1,0) \) lub \( {\bf n }=(0,1) \).
Z kolei nasz funkcjonał \( b(p ,{\bf v }) \) należy zmodyfikować dla metody DG w następujący sposób
\( b_{DG}(p, {\bf v})= - \sum_{K\in \Omega_h } \int_K \left ( p \frac{\partial v_1}{\partial x}+ p \frac{\partial v_2}{\partial y } \right) \)
\( + \sum_{F \in F_h } \int_F\left( [v_1] n_1 p+ [v_2] n_2 p \right) \)
Znaczenie symboli w całkach jest podobne jak w przypadku funkcjonału \( a({\bf u },{\bf v }) \).
Nasz stabilizowany problem Stokes'a wygląda teraz następująco:
\( a_{DG}({\bf u},{\bf v}) + b_{DG}(p,{\bf v }) = f({\bf v }) \quad \forall {\bf v} \)
\( b_{DG}(q,{\bf u } ) + {\color{red} {S_{DG}(p,q) }}= 0 \quad \forall q \)
Zgodnie z książką Daniele Antoni Di Pietro, Alexandre Ern, Mathematical Aspects of Discontinuous Galerkin Methods, w celu ustabilizowania równań Stokesa należy dodać tutaj człon
\( {\color{red}{S_{DG}(p,q)}} = \sum_{F\in F_h^{internal} } h_F \int_F [p][q] \)
w której \( F_h^{internal} \) oznacza wszystkie krawędzie wewnętrzne, tak więc sumujemy tutaj całki policzone po wszystkich krawędziach znajdujących się wewnątrz patcha elementów (omijamy krawędzie położone na brzegu). Podobnie jak poprzednio, \( h_F \) oznacza długość krawędzi.
W celu rozwiązania równań Stokesa stabilizowanych za pomocą metody DG, generujemy w podobny sposób układ równań, uwzględniając powyższe modyfikacje.
Możliwe jest również rozważanie różnych modyfikacji.
Założymy następujące fakty o funkcjach aproksymujących i testujących. Załóżmy że funkcje testujące będą kawałkami ciągłe, w sensie że nad poszczególnymi przedziałami wektora węzłów w analizie izogeometrycznej, lub nad poszczególnymi elementami skończonymi w sensie klasycznej metody elementów skończonych będą ciągłe, ale pomiędzy elemetami lub pomiędzy przedziałami wektora węzłów będą nieciągłe. Mamy więc "połamane" wielomiany testujące. Załóżmy też że funkcje aprokszmaczjne są ciągłe. Wówczas wszystkie człony zawierające skoki i średnie z funkcji aproksymujące znikną, w szczególności dostaniemy zmodyfikowane wyrażenia
\( a_{DG}({\bf u }, {\bf v }) = \sum_{K \in \Omega_h } \int_K \left( \frac{\partial u_1 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_1 }{\partial y } \frac{\partial v_1 }{\partial y } + \frac{\partial u_2 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_2 }{\partial y } \frac{\partial v_1 }{\partial y } \right) \)
\( -\sum_{F\in F_h} \int_F \left( \{\frac{\partial u_1}{\partial n_1} \} [\frac{\partial v_1}{\partial x}] + \{ \frac{\partial u_1}{\partial n_2} \} [\frac{\partial v_1}{\partial y } ] + \{\frac{\partial u_2}{\partial n_1} \} [\frac{\partial v_2}{\partial x}] + \{ \frac{\partial u_2}{\partial n_2} \} [\frac{\partial v_2}{\partial y } ]\right) \)
\( b_{DG}(p, {\bf v})= - \sum_{K\in \Omega_h } \int_K \left ( p \frac{\partial v_1}{\partial x}+ p \frac{\partial v_2}{\partial y } \right) \)
\( + \sum_{F \in F_h } \int_F\left( [v_1] n_1 p+ [v_2] n_2 p \right) \)
\( {\color{red}{S_{DG}(p,q)}} = 0 \)
Założymy następujące fakty o funkcjach aproksymujących i testujących. Załóżmy że funkcje aproksymujące będą kawałkami ciągłe, w sensie że nad poszczególnymi przedziałami wektora węzłów w analizie izogeometrycznej, lub nad poszczególnymi elementami skończonymi w sensie klasycznej metody elementów skończonych będą ciągłe, ale pomiędzy elemetami lub pomiędzy przedziałami wektora węzłów będą nieciągłe. Mamy więc "połamane" wielomiany aproksymacyjne. Załóżmy też że funkcje testujące są ciągłe. Wówczas mamy
\( a_{DG}({\bf u }, {\bf v }) = \sum_{K \in \Omega_h } \int_K \left( \frac{\partial u_1 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_1 }{\partial y } \frac{\partial v_1 }{\partial y } + \frac{\partial u_2 }{\partial x } \frac{\partial v_1}{\partial x } + \frac{\partial u_2 }{\partial y } \frac{\partial v_1 }{\partial y } \right) \)
\( - \sum_{F\in F_h} \int_F \left( [\frac{\partial u_1}{\partial x} ] \{ \frac{\partial v_1}{\partial n_1} \} + [\frac{\partial u_1}{\partial y} ] \{\frac{\partial v_1}{\partial n_2 } \} + [\frac{\partial u_2}{\partial x} ] \{ \frac{\partial v_2}{\partial n_1} \} + [\frac{\partial u_2}{\partial y} ] \{\frac{\partial v_2}{\partial n_2 } \} \right) \)
\( b_{DG}(p, {\bf v})= - \sum_{K\in \Omega_h } \int_K \left ( p \frac{\partial v_1}{\partial x}+ p \frac{\partial v_2}{\partial y } \right) \)
\( {\color{red}{S_{DG}(p,q)}} = 0 \)
Jeśli do sformułowania DG wsadzimy ciągłe przestrzenie aproksymacyjne i testujące, dostaniemy wówczas takie samo oryginalne sformułowanie Stokesa bez stabilizacji. Różne warianty przestrzeni aproksymacyjnych i testujących mogą być pomocne podczas testowania różnych form stabilizacji oraz do samego debugowania kodu generującego macierz metody DG.